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ABSTRACT 

We investigate the structure of the dark matter halo formed in three different 
cold dark matter scenarios. We performed N-body simulations of formation of 
13 cluster-sized halos. In all runs, density cusps proportional to r~^'^ developed 
at the center. This result was independent of the cosmological models we 
simulated. We could not reproduce the cusp shallower than r~^-^, which was 
obtained in some of previous studies. We also found that in all runs the density 
structure evolves in a self-similar way, even in Q 1 universes. These results 
show that the formation of structural form is a process decoupled from a 
background cosmology. 

Subject headings: cosmology:theory — dark matter — galaxies: clusters: general 
— methods: N-body simulations 



1. Introduction 

The structure of dark matter halos formed through disspationless hierarchical clustering 
from cosmological initial setting has been explored by many researchers since the "finding" 
of the universal profile by Navarro, Frenk and White (1996, 1997, hereafter NFW). NEW 
performed A^-body simulations of the halo formation and found that the density profiles of 
dark matter halos were expressed well by a simple formula. 



(r/ro)(l + r/ro)2 
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where po is a characteristic density and Tq is a scale radius. They argued that the profile of 
dark matter halos have the same shape, independent of the halo mass, the initial density 
fluctuation spectrum or the value of the cosmological parameters. 

Although the NFW results have been confirmed regarding the logarithmic slope 
shallower than isothermal near the center and —3 in the outskirt, disagreements concerning 
the slope of the central region were reported by subsequent studies in which higher-resolution 
simulations were performed. The disagreements are summarized into the following two 
kinds. 

One is that the slope at the center is steeper than that in the NFW result. Fukushige 
and Makino (1997, hereafter FM97) performed a simulation with 768k particles, while 
previous studies employed ~ 20k, and found that the galaxy-sized halo in their simulation 
has a cusp steeper than p oc r^^. Moore ct al. (1998, 1999 hereafter M99) and Ghigna ct 
al. (2000) performed simulations with up to 4M particles and obtained the results that the 
profile has a cusp proportional to r~^^^ both in galaxy-sized and cluster-sized halos. They 
proposed the modified universal profile (hereafter, the M99 profile), 

^ (r/r,)i-5[l + (r/r,)i-5]- 

Klypin et al. (2001) also obtained the results the the slope at the center can be approximated 
by r~^'^, though they argued that the NFW fit is still good up to their resolution limit. 

The second is that the slope at the center is not universal. Jing and Suto (2000) 
performed a series of A^-body simulations and concluded that the power of the cusp 
depends on mass, in contradiction to the claims by earlier studies. It varies from —1.5 for 
galaxy mass halo to —1.1 for cluster mass halo. Also, many analytical studies argued that 
the halo profile should depend on the power spectrum of initial density fluctuation. For 
example, Hoffman and Shaham (1985) and Syer and White (1998) predicted the slope, 
3(n -|- 3)/(n -|- 4) and 3(n -|- 3) / (n -|- 5) , respectively, where n is the effective power- law index 
of the power spectrum. 

In the first paper of this series (Fukushige and Makino 2001, hereafter Paper I), we 
performed simulations of 12 halos with the mass ranging from 6.6 x 10^^ M© to 8.0 x lO^^M©. 
We found that, in all runs, the halos have density cusps proportional to r~^"^ developed 
at the center. This result means that the density structure is universal in the sense that 
it is independent of the halo mass. We also found that the density structure evolves in a 
self-similar way as the central cusp grows outward with keeping the density near the center 
unchanged. 

This paper is a follow up of Paper I. In this paper, we focus on the dependence 
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of the halo profile on the cosmological model. We simulate halos in three conventional 
cold dark matter models: Standard, Lambda and Open CDM models, while simulations 
in Paper I was performed only in the standard CDM model. We performed A^-body 
simulations of formation of 13 cluster-sized dark matter halos using Barnes-Hut treecode 
and a special-purpose computer GRAPE-5 (Kawai et al. 2000). This is the first study 
that investigates the dependence on the different CDM models with the mass resolution 
and number of samples enough to discuss the slope of the cusp, though several simulations 
were performed in the same context (NFW, Thomas et al. 1998, Huss, Jain, and Steinmetz 
1999). 

The structure of this paper is as follows. In section 2, we describe the model of our 
A^-body simulation. In section 3, we present the results of simulation. Section 4 is for 
conclusion and discussion. 

2. Simulation Method 

We simulated the formation of the dark matter halos using "re-simulation" method, 
which has been the standard method for the simulation of halo formation since NFW. In 
this method, we first perform large scale cosmological simulations in order to identify the 
halo candidates. We trace back to the initial condition of the large scale simulation, and 
express the halo candidate with larger number of particles by adding a shorter wavelength 
perturbation. Then, we resimulate the halo candidates. 

We used three cosmological models listed in Table Standard, Open and Lambda 
Cold Dark Matter models (SCDM, OCDM, and LCDM). Here, fio is the density parameter, 
Ao is the dimensionless cosmological constant, Hq = 100/ikm ■ s ■ Mpc~^ at the present 
epoch, and Zi is the initial redshift. The amplitudes of the power spectrum in CDM models 
were normalized using the top-hat filtered mass variance at 8h~^ Mpc according to the 
cluster abundance (Kitayama & Suto 1997). 

The large scale cosmological simulations were performed with 1.1 x 10^ particles in a 
sphere of lOOMpc radius. The procedure for setting the initial condition were the same as 
those used in Fukushige and Suto (2001). We regard a spherical overdensity region around 
a local potential minimum within as a candidate halo. We define the radius such that 
the spherical overdensity inside is 178f2g^ times the critical density for SCDM and OCDM 
model, and 178^^-^ times for LCDM model (Eke, Cole, Frenk 1996). 

We selected 13 halos from the candidate halos catalog, which are summarized in Table 
0. We first selected the most massive halo in each CDM model, and then selected the halos 
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at random from halo candidates containing more than 1000 particles. We expressed a region 
within 5rv from the center of the halo at z = in the cosmological simulation with larger 
number of particles. We placed particles whose mass is the same as that in the cosmological 
simulation in a sphere of ~ 50Mpc radius surrounding the high resolution region in order 
to express the external tidal field. The total number of particles, A^, is listed in Table |[ 
As a result, the mass ratio between the high resolution particles and surrounding particles 
becomes rather large, 256 ~ 2048. In order to prevent the contamination of heavy particles 
into the halo, we set the boundary to be a rather large value, br^. No heavy particle entered 
within throughout all simulations. 

We integrated the system directly in the physical coordinates for both the cosmological 
and halo simulations (as in FM97, Paper I). We used a leap-fiog integrator with shared and 
constant timestep. The step size for the cosmological simulation is At/tn = 1/1024 and 
that for the halo simulation is listed in Table ^. The code for the time integration is the 
same as that in Fukushige and Suto (2001). We used the usual Plummer softening. The 
gravitational softening e is constant in the physical coordinates and the length is 5kpc for 
the cosmological simulation and Ikpc for all halo simulations. 

The force calculation is done with the Barnes- Hut tree code (Barnes & Hut 1986, 
Makino 1991) on GRAPE-5 (Kawai et al. 2000), a special-purpose computer designed to 
accelerate A^-body simulations. For most simulations, we used the GRAPE-5 system at the 
Astronomical Data Analysis Center of the National Astronomical Observatory, Japan. We 
used the opening parameter 9 = 0.4 for the cosmological simulation and 6 = 0.5 for the 
halo simulation. The simulations presented below required, for example in Run L6, ~ 80 
sees per timestep, and thus one run (16k timesteps) was completed in 370 CPU hours with 
a GRAPE-5 board connected to a host workstation with Alpha 21264 CPU (833MHz). 



Table 1: Cosmological Models 



Model 


^^0 


Ao 


h 


0-8 






SCDM 


1.0 


0.0 


0.5 


0.6 


6.9 X 10-8 


24.0 


LCDM 


0.3 


0.7 


0.7 


1.0 


1.4 X 10"^ 


32.3 


OCDM 


0.3 


0.0 


0.7 


1.0 


1.4 X 10-^ 


32.3 



3. Results 
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3.1. Snapshots 

Figure ^ shows the particle distribution for aU runs at 2; = 0. The length of the side 
for each panel is equal to 2r^. For these plots, we shifted the origin of coordinates to the 
position of the potential minimum. In Table ^, we summarized the radius r^, the mass M^, 
and the number of particles A^v within at z = 0. 

3.2. Density Profile 

Figure ^ shows the density profiles for all runs at z = 0. For Run 02, we plot the 
density profile at 2; = 0.026 because the merging process occurs just near the center of halo 
at 2; = 0. The position of the center of the halo was determined using the potential minimum 
and the density is averaged over each spherical shell whose width is log]^Q(Ar) = 0.0125. 
For the illustrative purpose, the densities are shifted vertically. 

In this figure, we plot the densities in the thick lines if the two criteria introduced in 
Paper I; (1) tj.ei{r)/t > 3 and (2) t^y{r)/At > 40, are satisfied, where tj.ei{r) is the local 
two-body relaxation time and tdy{r) is the local dynamical time. This means that the 
densities in the central region plotted in the thin lines are influenced by the numerical 
artifacts (mainly two-body relaxation), and they are not reliable. In the following discussion, 
we only use the densities plotted in the thick lines. 

In all runs we can see the central density cusps approximately proportional to r~^'^. In 
other words, the power of the cusp is —1.5 and is independent of cosmological models we 
simulated. The shallowing of the power-law index of the inner cusp observed in the LCDM 
runs by Jing and Suto (2000) was not reproduced in our LCDM runs. 

Moreover, the density profiles are in good agreement to the profile given by equation 
(|) (the M99 profile) in all runs. We set the scale radii tq as 0.8, 0.5, 0.4, 0.5 Mpc for S1..S4 
runs, 0.6, 0.5, 0.5, 0.35, 0.3, 0.3 Mpc for LI. .16 runs, and 0.35, 0.35, 0.25 Mpc for 01. .03 
runs. The agreement is very good for the inner region. In the outer region the agreement is 
not very good simply because the outer profile shows large fluctuations caused by individual 
infalling halos. The degree of the agreement is better in LCDM and OCDM model than in 
SCDM model. This is because the halo in LCDM and OCDM model is typically formed 
earlier and it is dynamically quiet around z ~ 0. 

Figure |^ shows the scale densities of the profile po [equation (j^)] and the concentration 
parameter c = r^/ro as a function of total mass. These values in S/LCDM models are 
consistent to those obtained by M99 and Jing and Suto (2000). We can see a tendency that 
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that cluster-size halos in OCDM models are more compact than that of their S/LCDM 
counterparts of the same mass. 



3.3. Self-Similarity 

In the SCDM simulations of Paper I, we found that the density profile evolves in a 
self-similar way as the central cusp grows outward keeping the density near the center 
unchanged (illustrated in Figure 14 of Paper I). However, it is not clear whether the density 
profile evolves in the same way in the LCDM and OCDM model, because neither of models 
retains the scale-free nature of the SCDM model. 

If the evolution is self-similar, we can write the density as 

p(r,t) = p^{M)p,{r,) (3) 
n = r/rt(M) (4) 
In Paper I, we obtained the self-similar variables, and r^, as 

n 

( M \^ , , 

'■•(^'^) = (lO^j 

assuming that the halo having a r" cusp grows outward in a self-similar way keeping the 
density of central cusp region constant and the fraction of mass in the cusp to total mass, 
M, is constant. Here n is the power of the cusp. 

In Figure § we plot the scaled density function of r,,. We plot three profiles 

at different value of the redshifts z for six runs. We set n to be —1.5 and use as the 
total mass. In this figure, we can see that the density profiles of the same halo at different 
times show good agreement to each other, which means that the density structure evolves 
self-similarly. Therefore, Figure ^ demonstrates that our assumption of the self-similarity is 
justified not only in the SCDM model, but also in LCDM and OCDM models. 

Finally, we consider whether halos of different masses are on the self-similar evolution 
track or not. If a halo evolves self-similarly, scaled scale densities p*o = Po/Pf should be 
constant during the evolution. Using equation (|), the scale density po is proportional to 

. In Figure ^, the dashed line indicates po oc . If the halo evolve self-similarly, the 
evolution track is along the dashed line in the — M plot. We can see a tendency that at 
least for 0/LCDM models, halos do not distribute along the self-similar evolution track, in 
other words, small halos cannot grow to large halos even if we wait longer. 
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4. Conclusion and Discussion 

We performed A^-body simulations of dark matter halo formation in three CDM 
models: Standard, Lambda, and Open CDM models. We simulated 13 halos whose mass 
range is 1.8 x lO^^M© to 8.6 x 10^^ M©. We used a widely adapted "re-simulation" method 
to set up initial conditions of halos and include the external tidal field. 

Our main conclusions are: 

(1) We found that, in all runs, the final halos have density cusps proportional to r~^'^, and 
the profiles show good agreement with the M99 profile, regardless of the cosmological 
models. 

(2) In all runs, the density profile evolves self-similarly. This is also independent of the 
cosmological models we simulated. 

There are some claims that the innermost slope should converge not to r^^-^ but to a 
shallower one (e.g. Taylor, Navarro 2001). Indeed, if we do not pay attention to numerical 
artifacts, we may see in our simulation results (Figure ^ that the innermost slopes of all 
runs become shallower than —1.5. However, Figure ^ shows that the inner region where the 
slope are shallower than r~^'^ is not reliable because in this region the accuracy criteria are 
not satisfied. The numerical artifact which makes the cusp shallower is mainly the two-body 
relaxation effect in this region (see Paper I). Therefore, we emphasize that discussions based 
on simulation results without careful analysis of the influence by the two-body relaxation 
effects are misleading. 

In our LCDM run, we could not reproduce the shallowing of the power-law index of 
the inner cusp observed in the LCDM runs by Jing and Suto (2000). This difference could 
be due to the smoothing by two-body relaxation in their cluster-sized halos. In this paper, 
we show that the density profile within ~ 0.01r2oo smoothed by the two-body relaxation. 
The density at 0.01r2oo and the mass resolution in their cluster-sized halo are similar to 
those in ours. The density profile in their simulations within 0.01r2oo; at which the profile 
begins to depart from r~^'^ inward, could be affected by the two-body relaxation. If their 
simulations had been performed with higher mass resolution, the slope might have approach 
to —1.5. The tendency can be already seen in their galaxy-sized halo. All halos have the 
cusp proportional to r~^'^ in their galaxy-sized halos whose central densities at 0.01r2oo are 
almost similar to those of their cluster-sized run and mass resolutions are about 100 times 
higher. 

Our results show that the NFW's claim concerning to the universality is certainly 
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valid. On the other hand, the density profile obtained are not in agreement with the NFW 
profile at the central region. Although it is important to find the convergence slope by 
further simulations, a simple explanation would be required for our final understanding of 
the universal profile. At present, we don't fully understand why the profile is universal 
and/or why the power of central cusp is —1.5, which we will address it in the future study. 
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hardware and software environment of the GRAPE-5 system, and to Yasushi Suto, Atsushi 
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by H.M.P.Couchman, P.A.Thomas, and F.R.Pearce. Most of numerical computations were 
carried out on the GRAPE system at ADAC (the Astronomical Data Analysis Center) of 
the National Astronomical Observatory, Japan. This research was partially supported by 
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no. JSPS-RFTP 97P01102. 
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Fig. 1. — Snapshots from all Runs ai z — 0. The length of the side for each panel is equal 
to 2ry. 
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Fig. 2. — Density profiles of tlie lialos for all Runs at ^ ~ 0. Only the densities plotted in 
the thick lines satisfy the accuracy criteria in section 3.2. Those plotted in the thin lines 
are influenced by numerical artifacts. The unit of density is Mq/pc^. The labels on the left 
of the profiles indicate the run name. The profiles except for Runs SI are vertically shifted 
downward by 1, 2, 3, 5, 6, 7 dex for Runs S2, S3, S4, 01, 02, 03, and by 1...6 for Runs 
L1...L6, respectively. The arrows indicate r^. The thin dashed lines indicate the densities 
proportional to r^^'^. The thin solid curves indicate the density profile given by equation 
(I) (the M99 profile). 
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Fig. 3. — Scale density po in the unit of Mq/pc^ and concentration parameter c = rv/ro as a 
function of total mass in solar mass. The star, square, and triangle symbols show those for 
SCDM, LCDM, and OCDM models, respectively. The dashed line indicates the self-similar 
evolutionary track discussed in section p73|. 
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Fig. 4. — Self-similar evolution of the density profile. The scaled densities arc plotted as 
a function of the scaled radius r^. The profile for Run S4, 01, 03, LI and L3 are vertically 
shifted downward by 1, 3, 4, 6, and 7 dex, respectively. The densities within 2ry at each 
redshift are plotted. 
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Table 2: Simulation Models 



Model 


Run 




Tv (Mpc) 




m (lO^Mo) 


Ai(lOV) 


N 


SCDM 


SI 


8.67 X 10^^ 


2.62 


1676525 


5.17 


1.57 


5074432 




S2 


5.46 X 10^^ 


2.21 


1056312 


5.17 


1.57 


3523844 




S3 


3.68 X lO^^ 


1.97 


1421930 


2.58 


1.57 


3478480 




S4 


3.58 X lO^^ 


1.91 


1383674 


2.58 


1.57 


4104120 


LCDM 


LI 


7.83 X lO^^ 


2.34 


1288779 


6.08 


0.82 


3624848 




L2 


5.32 X 10^^ 


2.08 


875058 


6.08 


1.63 


4360512 




L3 


3.97 X 10^4 


1.85 


1306187 


3.04 


0.82 


3066944 




L4 


2.17 X 10^^ 


1.52 


1425526 


1.52 


0.82 


3536640 




L5 


2.15 X 10^^ 


1.52 


707569 


3.04 


1.63 


2058140 




L6 


1.83 X 10^^ 


1.43 


1809105 


1.01 


0.82 


5458688 


OCDM 


01 


8.58 X 10^^ 


2.34 


1411523 


6.08 


0.68 


3711232 




02 


4.27 X 10^4 


1.86 


702022 


6.08 


1.37 


1748480 




03 


2.18 X 10^^ 


1.47 


717056 


3.04 


1.37 


2148736 



